home *** CD-ROM | disk | FTP | other *** search
/ Languguage OS 2 / Languguage OS II Version 10-94 (Knowledge Media)(1994).ISO / gnu / gmp-132.lha / gmp-1.3.2 / mpn_div.c < prev    next >
C/C++ Source or Header  |  1993-05-03  |  7KB  |  322 lines
  1. /* mpn_div -- Divide natural numbers, producing both remainder and
  2.    quotient.
  3.  
  4. Copyright (C) 1991 Free Software Foundation, Inc.
  5.  
  6. This file is part of the GNU MP Library.
  7.  
  8. The GNU MP Library is free software; you can redistribute it and/or modify
  9. it under the terms of the GNU General Public License as published by
  10. the Free Software Foundation; either version 2, or (at your option)
  11. any later version.
  12.  
  13. The GNU MP Library is distributed in the hope that it will be useful,
  14. but WITHOUT ANY WARRANTY; without even the implied warranty of
  15. MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  16. GNU General Public License for more details.
  17.  
  18. You should have received a copy of the GNU General Public License
  19. along with the GNU MP Library; see the file COPYING.  If not, write to
  20. the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.  */
  21.  
  22. #include "gmp.h"
  23. #include "gmp-impl.h"
  24. #include "longlong.h"
  25.  
  26. /* Divide num (NUM_PTR/NUM_SIZE) by den (DEN_PTR/DEN_SIZE) and write
  27.    the quotient at QUOT_PTR and the remainder at NUM_PTR.
  28.  
  29.    Return 0 or 1, depending on if the quotient size is (NSIZE - DSIZE)
  30.    or (NSIZE - DSIZE + 1).
  31.  
  32.    Argument constraints:
  33.    1. The most significant bit of d must be set.
  34.    2. QUOT_PTR != DEN_PTR and QUOT_PTR != NUM_PTR, i.e. the quotient storage
  35.       area must be distinct from either input operands.
  36.  
  37.    The exact sizes of the quotient and remainder must be determined
  38.    by the caller, in spite of the return value.  The return value just
  39.    informs the caller about if the highest digit is written or not, and
  40.    it may very well be 0.  */
  41.  
  42. /* THIS WILL BE IMPROVED SOON.  MORE COMMENTS AND FASTER CODE.  */
  43.  
  44. mp_size
  45. #ifdef __STDC__
  46. mpn_div (mp_ptr quot_ptr,
  47.      mp_ptr num_ptr, mp_size num_size,
  48.      mp_srcptr den_ptr, mp_size den_size)
  49. #else
  50. mpn_div (quot_ptr, num_ptr, num_size, den_ptr, den_size)
  51.      mp_ptr quot_ptr;
  52.      mp_ptr num_ptr;
  53.      mp_size num_size;
  54.      mp_srcptr den_ptr;
  55.      mp_size den_size;
  56. #endif
  57. {
  58.   mp_size q_is_long = 0;
  59.  
  60.   switch (den_size)
  61.     {
  62.     case 0:
  63.       /* We are asked to divide by zero, so go ahead and do it!
  64.      (To make the compiler not remove this statement, assign NUM_SIZE
  65.      and fall through.)  */
  66.       num_size = 1 / den_size;
  67.  
  68.     case 1:
  69.       {
  70.     mp_size i;
  71.     mp_limb n1, n0;
  72.     mp_limb d;
  73.  
  74.     d = den_ptr[0];
  75.     i = num_size - 1;
  76.     n1 = num_ptr[i];
  77.     i--;
  78.  
  79.     if (n1 >= d)
  80.       {
  81.         q_is_long = 1;
  82.         n1 = 0;
  83.         i++;
  84.       }
  85.  
  86.     for (; i >= 0; i--)
  87.       {
  88.         n0 = num_ptr[i];
  89.         udiv_qrnnd (quot_ptr[i], n1, n1, n0, d);
  90.       }
  91.  
  92.     num_ptr[0] = n1;
  93.       }
  94.       break;
  95.  
  96.     case 2:
  97.       {
  98.     mp_size i;
  99.     mp_limb n0, n1, n2;
  100.     mp_limb d0, d1;
  101.  
  102.     num_ptr += num_size - 2;
  103.     d0 = den_ptr[1];
  104.     d1 = den_ptr[0];
  105.     n0 = num_ptr[1];
  106.     n1 = num_ptr[0];
  107.  
  108.     if (n0 >= d0)
  109.       {
  110.         q_is_long = 1;
  111.         n1 = n0;
  112.         n0 = 0;
  113.         num_ptr++;
  114.         num_size++;
  115.       }
  116.  
  117.     for (i = num_size - den_size - 1; i >= 0; i--)
  118.       {
  119.         mp_limb q;
  120.         mp_limb r;
  121.  
  122.         num_ptr--;
  123.         if (n0 == d0)
  124.           {
  125.         /* Q should be either 111..111 or 111..110.  Need special
  126.            treatment of this rare case as normal division would
  127.            give overflow.  */
  128.         q = ~0;
  129.  
  130.         r = n1 + d0;
  131.         if (r < d0)    /* Carry in the addition? */
  132.           {
  133.             n2 = num_ptr[0];
  134.  
  135.             add_ssaaaa (n0, n1, r - d1, n2, 0, d1);
  136.             quot_ptr[i] = q;
  137.             continue;
  138.           }
  139.         n0 = d1 - (d1 != 0);
  140.         n1 = -d1;
  141.           }
  142.         else
  143.           {
  144.         udiv_qrnnd (q, r, n0, n1, d0);
  145.         umul_ppmm (n0, n1, d1, q);
  146.           }
  147.  
  148.         n2 = num_ptr[0];
  149.       q_test:
  150.         if (n0 > r || (n0 == r && n1 > n2))
  151.           {
  152.         /* The estimated Q was too large.  */
  153.         q--;
  154.  
  155.         sub_ddmmss (n0, n1, n0, n1, 0, d1);
  156.         r += d0;
  157.         if (r >= d0)    /* If not carry, test q again.  */
  158.           goto q_test;
  159.           }
  160.  
  161.         quot_ptr[i] = q;
  162.         sub_ddmmss (n0, n1, r, n2, n0, n1);
  163.       }
  164.     num_ptr[1] = n0;
  165.     num_ptr[0] = n1;
  166.       }
  167.       break;
  168.  
  169.     default:
  170.       {
  171.     mp_size i;
  172.     mp_limb d0 = den_ptr[den_size - 1];
  173.     mp_limb d1 = den_ptr[den_size - 2];
  174.     mp_limb n0 = num_ptr[num_size - 1];
  175.     int ugly_hack_flag = 0;
  176.  
  177.     if (n0 >= d0)
  178.       {
  179.  
  180.         /* There's a problem with this case, which shows up later in the
  181.            code.  q becomes 1 (and sometimes 0) the first time when
  182.            we've been here, and w_cy == 0 after the main do-loops below.
  183.            But c = num_ptr[j] reads rubbish outside the num_ptr vector!
  184.            Maybe I can solve this cleanly when I fix the early-end
  185.            optimization here in the default case.  For now, I change the
  186.            add_back entering condition, to kludge.  Leaving the stray
  187.            memref behind!
  188.  
  189.            HACK: Added ugly_hack_flag to make it work.  */
  190.  
  191.         q_is_long = 1;
  192.         n0 = 0;
  193.         num_size++;
  194.         ugly_hack_flag = 1;
  195.       }
  196.  
  197.     num_ptr += num_size;
  198.     den_ptr += den_size;
  199.     for (i = num_size - den_size - 1; i >= 0; i--)
  200.       {
  201.         mp_limb q;
  202.         mp_limb n1;
  203.         mp_limb w_cy;
  204.         mp_limb d, c;
  205.         mp_size j;
  206.  
  207.         num_ptr--;
  208.         if (n0 == d0)
  209.           /* This might over-estimate q, but it's probably not worth
  210.          the extra code here to find out.  */
  211.           q = ~0;
  212.         else
  213.           {
  214.         mp_limb r;
  215.  
  216.         udiv_qrnnd (q, r, n0, num_ptr[-1], d0);
  217.         umul_ppmm (n1, n0, d1, q);
  218.  
  219.         while (n1 > r || (n1 == r && n0 > num_ptr[-2]))
  220.           {
  221.             q--;
  222.             r += d0;
  223.             if (r < d0)    /* I.e. "carry in previous addition?"  */
  224.               break;
  225.             n1 -= n0 < d1;
  226.             n0 -= d1;
  227.           }
  228.           }
  229.  
  230.         w_cy = 0;
  231.         j = -den_size;
  232.         do
  233.           {
  234.         d = den_ptr[j];
  235.         c = num_ptr[j];
  236.         umul_ppmm (n1, n0, d, q);
  237.         n0 += w_cy;
  238.         w_cy = (n0 < w_cy) + n1;
  239.         n0 = c - n0;
  240.         num_ptr[j] = n0;
  241.         if (n0 > c)
  242.           goto cy_loop;
  243.           ncy_loop:
  244.         j++;
  245.           }
  246.         while  (j < 0);
  247.  
  248.         if (ugly_hack_flag)
  249.           {
  250.         c = 0;
  251.         ugly_hack_flag = 0;
  252.           }
  253.         else
  254.           c = num_ptr[j];
  255.         if (c >= w_cy)
  256.           goto store_q;
  257.         goto add_back;
  258.  
  259.         do
  260.           {
  261.         d = den_ptr[j];
  262.         c = num_ptr[j];
  263.         umul_ppmm (n1, n0, d, q);
  264.         n0 += w_cy;
  265.         w_cy = (n0 < w_cy) + n1;
  266.         n0 = c - n0 - 1;
  267.         num_ptr[j] = n0;
  268.         if (n0 < c)
  269.           goto ncy_loop;
  270.           cy_loop:
  271.         j++;
  272.           }
  273.         while  (j < 0);
  274.  
  275.         if (ugly_hack_flag)
  276.           {
  277.         c = 0;
  278.         ugly_hack_flag = 0;
  279.           }
  280.         else
  281.           c = num_ptr[j];
  282.         w_cy++;
  283.         if (c >= w_cy)
  284.           goto store_q;
  285.  
  286.       add_back:
  287.         j = -den_size;
  288.         do
  289.           {
  290.         d = den_ptr[j];
  291.         n0 = num_ptr[j] + d;
  292.         num_ptr[j] = n0;
  293.         if (n0 < d)
  294.           goto ab_cy_loop;
  295.           ab_ncy_loop:
  296.         j++;
  297.           }
  298.         while  (j < 0);
  299.         abort ();        /* We should always have a carry out! */
  300.  
  301.         do
  302.           {
  303.         d = den_ptr[j];
  304.         n0 = num_ptr[j] + d + 1;
  305.         num_ptr[j] = n0;
  306.         if (n0 > d)
  307.           goto ab_ncy_loop;
  308.           ab_cy_loop:
  309.         j++;
  310.           }
  311.         while  (j < 0);
  312.         q--;
  313.  
  314.       store_q:
  315.         quot_ptr[i] = q;
  316.       }
  317.       }
  318.     }
  319.  
  320.   return q_is_long;
  321. }
  322.